source('misc_functions/functions.R')
library(tidyverse)
library(plotly)
library(modelr)
library(mgcv)
The data set has been constructed using average Science scores by country from the Programme for International Student Assessment (PISA) 2006, along with GNI per capita (Purchasing Power Parity, 2005 dollars), Educational Index, Health Index, and Human Development Index from UN data. The key variables are as follows:
- Overall Science Score (average score for 15 year olds)
- Interest in science
- Support for scientific inquiry
- Income Index
- Health Index
- Education Index
- Human Development Index (composed of the Income index, Health Index, and Education Index)
The first thing to do is get the data in and do some initial inspections.
We can look at the individual relationships of covariates with overall science score.
dmelt = pisa %>%
select(-Evidence, -Explain, -Issues) %>%
gather(key=Variable,
value=Value,
-Overall, -Country)
ggplot(aes(x=Value,y=Overall), data=dmelt) +
geom_point(color='#ff5500',alpha=.75) +
geom_smooth(se=F, lwd=.5, color='#00aaff') +
geom_text(aes(label=Country), alpha=0, size=1,angle=30, hjust=-.2, # making transparent so only plotly will show the country
vjust=-.2) +
facet_wrap(~Variable, scales='free_x') +
labs(x='') +
theme_trueMinimal()

Single Predictor
Linear Fit
We will start with the simple situation of a single predictor. Let’s begin by using a typical linear regression to predict science scores by the Income index. We could use the standard R lm function, but I’ll leave that as an exercise for comparison. We can still do straightforward linear models with the gam function, and again it is important to note that the standard linear model can be seen as a special case of a GAM.
Family: gaussian
Link function: identity
Formula:
Overall ~ Income
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 204.32 35.37 5.777 4.32e-07 ***
Income 355.85 46.79 7.606 5.36e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.518 Deviance explained = 52.7%
GCV = 1504.5 Scale est. = 1448.8 n = 54
GAM
Fitting the model
Let’s now fit an actual generalized additive model using the same cubic spline as our basis. We again use the gam function as before for basic model fitting, but now we are using a function s within the formula to denote the smooth terms. Within that function we also specify the type of smooth, though a default is available. I chose bs = cr, denoting cubic regression splines.
mod_gam1 <- gam(Overall ~ s(Income, bs="cr"), data=pisa)
summary(mod_gam1)
Family: gaussian
Link function: identity
Formula:
Overall ~ s(Income, bs = "cr")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 470.444 4.082 115.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Income) 6.895 7.741 16.67 1.59e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.7 Deviance explained = 73.9%
GCV = 1053.7 Scale est. = 899.67 n = 54
In the summary, we first see the distribution assumed, as well as the link function used, in this case normal and identity, respectively, which to iterate, had we had no smoothing, would result in a standard linear model.
We have two general pieces of output:
- Parametric: any non-smooth term, interpreted as any coefficient in the same GLM model.
- Smooth terms: stuff we allow to ‘wiggle’ (or otherwise have a special relationship), interpreted visually via component plot.
In this case, the only parametric component is the intercept, but it’s good to remember that you are not bound to smooth every effect of interest, and indeed, as we will discuss in more detail later, part of the process may involve refitting the model with terms that were found to be linear for the most part anyway. The smooth component of our model regarding a country’s income and its relationship with overall science score suggests that it is statistically significant, but there are a couple of things in the model summary that would be unfamiliar. You may consult the primary document for details, but the main thing is that you can interpret it like any result you’d see in a regression table. The statistical significance is an approximation though, so again, let visualization be your guide.
Since this is in the end a linear model with a normal distribution for the target, we can use the R^2 value as we would with lm. Note that it already is adjusted for you. The deviance explained is the unadjusted R^2 in this case, but generalizes to other models that don’t assume the normal distribution. The scale estimate is our residual sums of squares from the standard regression setting, and the GCV is a better estimate of what that would be on held-out (i.e. new) data. The point is that you are seeing the usual output you have with other models, just with different labels.
Graphical Display
The primary way to interpret the effect of income is visually though. The mgcv package will provide a basic plot as follows:
plot(mod_gam1)

What you’re seeing is a component plot with the y/fitted values centered. Almost no one finds this intuitive in practice, nor is it usually what they’d report. Instead we can get predictions at key values of other covariates (e.g. held at their mean or reference category). In this case with only one predictor, we just get the predicted values on the original scale of the outcome. We can use the ggeffects package for this (or see my visibly package).
library(ggeffects)
plot_dat <- ggpredict(mod_gam1, terms = "Income")
ggplot(plot_dat, aes(x = x, y = predicted)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
geom_line(color = 'dodgerblue') +
labs(x = 'Income')

Model Comparison
We can compare models via AIC or the GCV coefficient, both of which focus on the predictive capabilities of the model on new data.
AIC(mod_lm, mod_gam1)
Likelihood ratio test (approximate).
anova(mod_lm, mod_gam1, test="Chisq")
Analysis of Deviance Table
Model 1: Overall ~ Income
Model 2: Overall ~ s(Income, bs = "cr")
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 52.000 75336
2 45.259 41479 6.7411 33857 2.778e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Multiple Predictors
Linear FIt
mod_lm2 <- gam(Overall ~ Income + Edu + Health, data=pisa)
summary(mod_lm2)
Family: gaussian
Link function: identity
Formula:
Overall ~ Income + Edu + Health
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 121.18 78.97 1.535 0.1314
Income 182.32 85.27 2.138 0.0376 *
Edu 234.11 54.78 4.274 9.06e-05 ***
Health 27.01 134.90 0.200 0.8421
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.616 Deviance explained = 63.9%
GCV = 1212.3 Scale est. = 1119 n = 52
GAM
mod_gam2 <- gam(Overall ~ s(Income) + s(Edu) + s(Health), data=pisa)
summary(mod_gam2)
Family: gaussian
Link function: identity
Formula:
Overall ~ s(Income) + s(Edu) + s(Health)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 471.154 2.772 170 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(Income) 7.593 8.415 8.826 1.33e-07 ***
s(Edu) 6.204 7.178 3.308 0.00733 **
s(Health) 1.000 1.000 2.736 0.10661
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.863 Deviance explained = 90.3%
GCV = 573.83 Scale est. = 399.5 n = 52
We can see the fit appears to be better. In addition, as the effective degrees of freedom is 1 for Health, it is essentially a linear fit, so if desired, we can leave it among the parametric terms.
# plot(mod_gam2)
library(patchwork) # devtools::install_github("thomasp85/patchwork")
g1 =
ggpredict(mod_gam2, terms = "Income") %>%
ggplot(aes(x = x, y = predicted)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
geom_line(color = 'dodgerblue') +
labs(x = 'Income')
g2 =
ggpredict(mod_gam2, terms = "Edu") %>%
ggplot(aes(x = x, y = predicted)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
geom_line(color = 'dodgerblue') +
labs(x = 'Edu')
g3 =
ggpredict(mod_gam2, terms = "Health") %>%
ggplot(aes(x = x, y = predicted)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
geom_line(color = 'dodgerblue') +
labs(x = 'Health')
(g2 + g3 + g1 + plot_layout(nrow = 2)) * theme_trueMinimal()

2D Smooths
mod_gam3 <- gam(Overall ~ Health + te(Income, Edu), data=pisa)
summary(mod_gam3)
Family: gaussian
Link function: identity
Formula:
Overall ~ Health + te(Income, Edu)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 574.5 101.4 5.666 1.42e-06 ***
Health -115.8 113.5 -1.020 0.314
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
te(Income,Edu) 10.31 12.41 9.084 2.43e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.804 Deviance explained = 84.8%
GCV = 747.52 Scale est. = 570.59 n = 52
# requires the visibly package, otherwise see the main document;
# devtools::install_github("m-clark/visibly")
visibly::plot_gam_3d(model = mod_gam3, main_var = Income, second_var = Edu, palette='bilbao', direction = 1)
LS0tCnRpdGxlOiAiQXBwbGljYXRpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawplZGl0b3Jfb3B0aW9uczogCiAgY2h1bmtfb3V0cHV0X3R5cGU6IGlubGluZQotLS0KCmBgYHtyIG1pc2NfZnVuY3Rpb25zfQpzb3VyY2UoJ21pc2NfZnVuY3Rpb25zL2Z1bmN0aW9ucy5SJykKYGBgCgpgYGB7ciBsb2FkX3BhY2thZ2VzLCBtZXNzYWdlPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShwbG90bHkpCmxpYnJhcnkobW9kZWxyKQpsaWJyYXJ5KG1nY3YpCmBgYAoKVGhlIGRhdGEgc2V0IGhhcyBiZWVuIGNvbnN0cnVjdGVkIHVzaW5nIGF2ZXJhZ2UgU2NpZW5jZSBzY29yZXMgYnkgY291bnRyeSBmcm9tIHRoZSBQcm9ncmFtbWUgZm9yIEludGVybmF0aW9uYWwgU3R1ZGVudCBBc3Nlc3NtZW50IChQSVNBKSAyMDA2LCBhbG9uZyB3aXRoIEdOSSBwZXIgY2FwaXRhIChQdXJjaGFzaW5nIFBvd2VyIFBhcml0eSwgMjAwNSBkb2xsYXJzKSwgRWR1Y2F0aW9uYWwgSW5kZXgsIEhlYWx0aCBJbmRleCwgYW5kIEh1bWFuIERldmVsb3BtZW50IEluZGV4IGZyb20gVU4gZGF0YS4gVGhlIGtleSB2YXJpYWJsZXMgYXJlIGFzIGZvbGxvd3M6CgoKLSBPdmVyYWxsIFNjaWVuY2UgU2NvcmUgKGF2ZXJhZ2Ugc2NvcmUgZm9yIDE1IHllYXIgb2xkcykKLSBJbnRlcmVzdCBpbiBzY2llbmNlCi0gU3VwcG9ydCBmb3Igc2NpZW50aWZpYyBpbnF1aXJ5Ci0gSW5jb21lIEluZGV4Ci0gSGVhbHRoIEluZGV4Ci0gRWR1Y2F0aW9uIEluZGV4Ci0gSHVtYW4gRGV2ZWxvcG1lbnQgSW5kZXggKGNvbXBvc2VkIG9mIHRoZSBJbmNvbWUgaW5kZXgsIEhlYWx0aCBJbmRleCwgYW5kIEVkdWNhdGlvbiBJbmRleCkKClRoZSBmaXJzdCB0aGluZyB0byBkbyBpcyBnZXQgdGhlIGRhdGEgaW4gYW5kIGRvIHNvbWUgaW5pdGlhbCBpbnNwZWN0aW9ucy4KCmBgYHtyIGluaXRpYWxfaW5zcGVjdGlvbl9vZl9waXNhLCBlY2hvPTF9CnBpc2EgPSByZWFkLmNzdignZGF0YS9waXNhc2NpMjAwNi5jc3YnKQoKcGlzYSAKCnBpc2EgJT4lIHNlbGVjdCgtQ291bnRyeSkKYGBgCgpXZSBjYW4gbG9vayBhdCB0aGUgaW5kaXZpZHVhbCByZWxhdGlvbnNoaXBzIG9mIGNvdmFyaWF0ZXMgd2l0aCBvdmVyYWxsIHNjaWVuY2Ugc2NvcmUuCgpgYGB7ciBiaXZhcmlhdGVfcmVsYXRpb25zaGlwcywgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRX0KZG1lbHQgPSBwaXNhICU+JSAKICBzZWxlY3QoLUV2aWRlbmNlLCAtRXhwbGFpbiwgLUlzc3VlcykgJT4lIAogIGdhdGhlcihrZXk9VmFyaWFibGUsIAogICAgICAgICB2YWx1ZT1WYWx1ZSwgCiAgICAgICAgIC1PdmVyYWxsLCAtQ291bnRyeSkKCmdncGxvdChhZXMoeD1WYWx1ZSx5PU92ZXJhbGwpLCBkYXRhPWRtZWx0KSArCiAgZ2VvbV9wb2ludChjb2xvcj0nI2ZmNTUwMCcsYWxwaGE9Ljc1KSArCiAgZ2VvbV9zbW9vdGgoc2U9RiwgbHdkPS41LCBjb2xvcj0nIzAwYWFmZicpICsKICBnZW9tX3RleHQoYWVzKGxhYmVsPUNvdW50cnkpLCBhbHBoYT0wLCBzaXplPTEsYW5nbGU9MzAsIGhqdXN0PS0uMiwgICAjIG1ha2luZyB0cmFuc3BhcmVudCBzbyBvbmx5IHBsb3RseSB3aWxsIHNob3cgdGhlIGNvdW50cnkKICAgICAgICAgICAgdmp1c3Q9LS4yKSArCiAgZmFjZXRfd3JhcCh+VmFyaWFibGUsIHNjYWxlcz0nZnJlZV94JykgKwogIGxhYnMoeD0nJykgKwogIHRoZW1lX3RydWVNaW5pbWFsKCkKYGBgCgoKIyMgU2luZ2xlIFByZWRpY3RvcgoKIyMjIExpbmVhciBGaXQKCldlIHdpbGwgc3RhcnQgd2l0aCB0aGUgc2ltcGxlIHNpdHVhdGlvbiBvZiBhIHNpbmdsZSBwcmVkaWN0b3IuIExldCdzIGJlZ2luIGJ5IHVzaW5nIGEgdHlwaWNhbCBsaW5lYXIgcmVncmVzc2lvbiB0byBwcmVkaWN0IHNjaWVuY2Ugc2NvcmVzIGJ5IHRoZSBJbmNvbWUgaW5kZXguICBXZSBjb3VsZCB1c2UgdGhlIHN0YW5kYXJkIFIgYGxtYCBmdW5jdGlvbiwgYnV0IEknbGwgbGVhdmUgdGhhdCBhcyBhbiBleGVyY2lzZSBmb3IgY29tcGFyaXNvbi4gIFdlIGNhbiBzdGlsbCBkbyBzdHJhaWdodGZvcndhcmQgbGluZWFyIG1vZGVscyB3aXRoIHRoZSBgZ2FtYCBmdW5jdGlvbiwgYW5kIGFnYWluIGl0IGlzIGltcG9ydGFudCB0byBub3RlIHRoYXQgdGhlIHN0YW5kYXJkIGxpbmVhciBtb2RlbCBjYW4gYmUgc2VlbiBhcyBhIHNwZWNpYWwgY2FzZSBvZiBhIEdBTS4KCmBgYHtyIG1vZF9sbSwgZWNobz0tNH0KbGlicmFyeShtZ2N2KQptb2RfbG0gPC0gZ2FtKE92ZXJhbGwgfiBJbmNvbWUsIGRhdGE9cGlzYSkKc3VtbWFyeShtb2RfbG0pCmBgYAoKIyMjIEdBTQoKIyMjIEZpdHRpbmcgdGhlIG1vZGVsCgpMZXQncyBub3cgZml0IGFuIGFjdHVhbCBnZW5lcmFsaXplZCBhZGRpdGl2ZSBtb2RlbCB1c2luZyB0aGUgc2FtZSBjdWJpYyBzcGxpbmUgYXMgb3VyIGJhc2lzLiBXZSBhZ2FpbiB1c2UgdGhlIGBnYW1gIGZ1bmN0aW9uIGFzIGJlZm9yZSBmb3IgYmFzaWMgbW9kZWwgZml0dGluZywgYnV0IG5vdyB3ZSBhcmUgdXNpbmcgYSBmdW5jdGlvbiBgc2Agd2l0aGluIHRoZSBmb3JtdWxhIHRvIGRlbm90ZSB0aGUgc21vb3RoIHRlcm1zLiAgV2l0aGluIHRoYXQgZnVuY3Rpb24gd2UgYWxzbyBzcGVjaWZ5IHRoZSB0eXBlIG9mIHNtb290aCwgdGhvdWdoIGEgZGVmYXVsdCBpcyBhdmFpbGFibGUuICBJIGNob3NlIGBicyA9IGNyYCwgZGVub3RpbmcgY3ViaWMgcmVncmVzc2lvbiBzcGxpbmVzLgoKYGBge3IgbW9kX2dhbTF9Cm1vZF9nYW0xIDwtIGdhbShPdmVyYWxsIH4gcyhJbmNvbWUsIGJzPSJjciIpLCBkYXRhPXBpc2EpCnN1bW1hcnkobW9kX2dhbTEpCmBgYAoKSW4gdGhlIHN1bW1hcnksIHdlIGZpcnN0IHNlZSB0aGUgZGlzdHJpYnV0aW9uIGFzc3VtZWQsIGFzIHdlbGwgYXMgdGhlIGxpbmsgZnVuY3Rpb24gdXNlZCwgaW4gdGhpcyBjYXNlIG5vcm1hbCBhbmQgaWRlbnRpdHksIHJlc3BlY3RpdmVseSwgd2hpY2ggdG8gaXRlcmF0ZSwgaGFkIHdlIGhhZCBubyBzbW9vdGhpbmcsIHdvdWxkIHJlc3VsdCBpbiBhIHN0YW5kYXJkIGxpbmVhciBtb2RlbC4KCldlIGhhdmUgdHdvIGdlbmVyYWwgcGllY2VzIG9mIG91dHB1dDoKCi0gKipQYXJhbWV0cmljKio6IGFueSBub24tc21vb3RoIHRlcm0sIGludGVycHJldGVkIGFzIGFueSBjb2VmZmljaWVudCBpbiB0aGUgc2FtZSBHTE0gbW9kZWwuCi0gKipTbW9vdGggdGVybXMqKjogc3R1ZmYgd2UgYWxsb3cgdG8gJ3dpZ2dsZScgKG9yIG90aGVyd2lzZSBoYXZlIGEgc3BlY2lhbCByZWxhdGlvbnNoaXApLCBpbnRlcnByZXRlZCB2aXN1YWxseSB2aWEgY29tcG9uZW50IHBsb3QuCgpJbiB0aGlzIGNhc2UsIHRoZSBvbmx5IHBhcmFtZXRyaWMgY29tcG9uZW50IGlzIHRoZSBpbnRlcmNlcHQsIGJ1dCBpdCdzIGdvb2QgdG8gcmVtZW1iZXIgdGhhdCB5b3UgYXJlIG5vdCBib3VuZCB0byBzbW9vdGggZXZlcnkgZWZmZWN0IG9mIGludGVyZXN0LCBhbmQgaW5kZWVkLCBhcyB3ZSB3aWxsIGRpc2N1c3MgaW4gbW9yZSBkZXRhaWwgbGF0ZXIsIHBhcnQgb2YgdGhlIHByb2Nlc3MgbWF5IGludm9sdmUgcmVmaXR0aW5nIHRoZSBtb2RlbCB3aXRoIHRlcm1zIHRoYXQgd2VyZSBmb3VuZCB0byBiZSBsaW5lYXIgZm9yIHRoZSBtb3N0IHBhcnQgYW55d2F5LiAgVGhlIHNtb290aCBjb21wb25lbnQgb2Ygb3VyIG1vZGVsIHJlZ2FyZGluZyBhIGNvdW50cnkncyBpbmNvbWUgYW5kIGl0cyByZWxhdGlvbnNoaXAgd2l0aCBvdmVyYWxsIHNjaWVuY2Ugc2NvcmUgc3VnZ2VzdHMgdGhhdCBpdCBpcyBzdGF0aXN0aWNhbGx5IHNpZ25pZmljYW50LCBidXQgdGhlcmUgYXJlIGEgY291cGxlIG9mIHRoaW5ncyBpbiB0aGUgbW9kZWwgc3VtbWFyeSB0aGF0IHdvdWxkIGJlIHVuZmFtaWxpYXIuICBZb3UgbWF5IGNvbnN1bHQgdGhlIHByaW1hcnkgZG9jdW1lbnQgZm9yIGRldGFpbHMsIGJ1dCB0aGUgbWFpbiB0aGluZyBpcyB0aGF0IHlvdSBjYW4gaW50ZXJwcmV0IGl0IGxpa2UgYW55IHJlc3VsdCB5b3UnZCBzZWUgaW4gYSByZWdyZXNzaW9uIHRhYmxlLiAgVGhlIHN0YXRpc3RpY2FsIHNpZ25pZmljYW5jZSBpcyBhbiBhcHByb3hpbWF0aW9uIHRob3VnaCwgc28gYWdhaW4sIGxldCB2aXN1YWxpemF0aW9uIGJlIHlvdXIgZ3VpZGUuCgpTaW5jZSB0aGlzIGlzIGluIHRoZSBlbmQgYSBsaW5lYXIgbW9kZWwgd2l0aCBhIG5vcm1hbCBkaXN0cmlidXRpb24gZm9yIHRoZSB0YXJnZXQsIHdlIGNhbiB1c2UgdGhlIFJeMiB2YWx1ZSBhcyB3ZSB3b3VsZCB3aXRoIGBsbWAuICBOb3RlIHRoYXQgaXQgYWxyZWFkeSBpcyBhZGp1c3RlZCBmb3IgeW91LiBUaGUgKmRldmlhbmNlIGV4cGxhaW5lZCogaXMgdGhlIHVuYWRqdXN0ZWQgUl4yIGluIHRoaXMgY2FzZSwgYnV0IGdlbmVyYWxpemVzIHRvIG90aGVyIG1vZGVscyB0aGF0IGRvbid0IGFzc3VtZSB0aGUgbm9ybWFsIGRpc3RyaWJ1dGlvbi4gVGhlICpzY2FsZSBlc3RpbWF0ZSogaXMgb3VyIHJlc2lkdWFsIHN1bXMgb2Ygc3F1YXJlcyBmcm9tIHRoZSBzdGFuZGFyZCByZWdyZXNzaW9uIHNldHRpbmcsIGFuZCB0aGUgKkdDViogaXMgYSBiZXR0ZXIgZXN0aW1hdGUgb2Ygd2hhdCB0aGF0IHdvdWxkIGJlIG9uIGhlbGQtb3V0IChpLmUuIG5ldykgZGF0YS4gIFRoZSBwb2ludCBpcyB0aGF0IHlvdSBhcmUgc2VlaW5nIHRoZSB1c3VhbCBvdXRwdXQgeW91IGhhdmUgd2l0aCBvdGhlciBtb2RlbHMsIGp1c3Qgd2l0aCBkaWZmZXJlbnQgbGFiZWxzLgoKIyMjIEdyYXBoaWNhbCBEaXNwbGF5CgpUaGUgcHJpbWFyeSB3YXkgdG8gaW50ZXJwcmV0IHRoZSBlZmZlY3Qgb2YgaW5jb21lIGlzIHZpc3VhbGx5IHRob3VnaC4gIFRoZSAqbWdjdiogcGFja2FnZSB3aWxsIHByb3ZpZGUgYSBiYXNpYyBwbG90IGFzIGZvbGxvd3M6CgpgYGB7ciBtZ2N2X3Bsb3R9CnBsb3QobW9kX2dhbTEpCmBgYAoKCldoYXQgeW91J3JlIHNlZWluZyBpcyBhICpjb21wb25lbnQgcGxvdCogd2l0aCB0aGUgeS9maXR0ZWQgdmFsdWVzIGNlbnRlcmVkLiAgQWxtb3N0IG5vIG9uZSBmaW5kcyB0aGlzIGludHVpdGl2ZSBpbiBwcmFjdGljZSwgbm9yIGlzIGl0IHVzdWFsbHkgd2hhdCB0aGV5J2QgcmVwb3J0LiAgSW5zdGVhZCB3ZSBjYW4gZ2V0IHByZWRpY3Rpb25zIGF0IGtleSB2YWx1ZXMgb2Ygb3RoZXIgY292YXJpYXRlcyAoZS5nLiBoZWxkIGF0IHRoZWlyIG1lYW4gb3IgcmVmZXJlbmNlIGNhdGVnb3J5KS4gIEluIHRoaXMgY2FzZSB3aXRoIG9ubHkgb25lIHByZWRpY3Rvciwgd2UganVzdCBnZXQgdGhlIHByZWRpY3RlZCB2YWx1ZXMgb24gdGhlIG9yaWdpbmFsIHNjYWxlIG9mIHRoZSBvdXRjb21lLiAgV2UgY2FuIHVzZSB0aGUgYGdnZWZmZWN0c2AgcGFja2FnZSBmb3IgdGhpcyAob3Igc2VlIG15IFt2aXNpYmx5IHBhY2thZ2VdKGh0dHBzOi8vbS1jbGFyay5naXRodWIuaW8vdmlzaWJseSkpLgoKYGBge3IgdmlzdWFsaXplX2luY29tZV9tYXJnaW5hbF9lZmZlY3R9CmxpYnJhcnkoZ2dlZmZlY3RzKQoKcGxvdF9kYXQgPC0gZ2dwcmVkaWN0KG1vZF9nYW0xLCB0ZXJtcyA9ICJJbmNvbWUiKQoKZ2dwbG90KHBsb3RfZGF0LCBhZXMoeCA9IHgsIHkgPSBwcmVkaWN0ZWQpKSArIAogIGdlb21fcmliYm9uKGFlcyh5bWluID0gY29uZi5sb3csIHltYXggPSBjb25mLmhpZ2gpLCBhbHBoYSA9IC4yNSkgKwogIGdlb21fbGluZShjb2xvciA9ICdkb2RnZXJibHVlJykgKyAKICBsYWJzKHggPSAnSW5jb21lJykKYGBgCgojIyMgTW9kZWwgQ29tcGFyaXNvbgoKV2UgY2FuIGNvbXBhcmUgbW9kZWxzIHZpYSBBSUMgb3IgdGhlIEdDViBjb2VmZmljaWVudCwgYm90aCBvZiB3aGljaCBmb2N1cyBvbiB0aGUgcHJlZGljdGl2ZSBjYXBhYmlsaXRpZXMgb2YgdGhlIG1vZGVsIG9uIG5ldyBkYXRhLgoKYGBge3IgbW9kZWxfY29tcGFyaXNvbn0KQUlDKG1vZF9sbSwgbW9kX2dhbTEpCmBgYAoKTGlrZWxpaG9vZCByYXRpbyB0ZXN0IChhcHByb3hpbWF0ZSkuCgpgYGB7ciBhbm92YV9nYW19CmFub3ZhKG1vZF9sbSwgbW9kX2dhbTEsIHRlc3Q9IkNoaXNxIikKYGBgCgoKIyMgTXVsdGlwbGUgUHJlZGljdG9ycwoKCiMjIyBMaW5lYXIgRkl0CgpgYGB7ciBtb2RfbG0yfQptb2RfbG0yIDwtIGdhbShPdmVyYWxsIH4gSW5jb21lICsgRWR1ICsgSGVhbHRoLCBkYXRhPXBpc2EpCnN1bW1hcnkobW9kX2xtMikKYGBgCgojIyMgR0FNCgpgYGB7ciBtb2RfZ2FtMn0KbW9kX2dhbTIgPC0gZ2FtKE92ZXJhbGwgfiBzKEluY29tZSkgKyBzKEVkdSkgKyBzKEhlYWx0aCksIGRhdGE9cGlzYSkKc3VtbWFyeShtb2RfZ2FtMikKYGBgCgpXZSBjYW4gc2VlIHRoZSBmaXQgYXBwZWFycyB0byBiZSBiZXR0ZXIuICBJbiBhZGRpdGlvbiwgYXMgdGhlIGVmZmVjdGl2ZSBkZWdyZWVzIG9mIGZyZWVkb20gaXMgMSBmb3IgSGVhbHRoLCBpdCBpcyBlc3NlbnRpYWxseSBhIGxpbmVhciBmaXQsIHNvIGlmIGRlc2lyZWQsIHdlIGNhbiBsZWF2ZSBpdCBhbW9uZyB0aGUgYHBhcmFtZXRyaWNgIHRlcm1zLgoKYGBge3IgbW9kX2dhbTJfcGxvdH0KIyBwbG90KG1vZF9nYW0yKQoKbGlicmFyeShwYXRjaHdvcmspICMgZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJ0aG9tYXNwODUvcGF0Y2h3b3JrIikKCmcxID0gCiAgZ2dwcmVkaWN0KG1vZF9nYW0yLCB0ZXJtcyA9ICJJbmNvbWUiKSAlPiUgCiAgZ2dwbG90KGFlcyh4ID0geCwgeSA9IHByZWRpY3RlZCkpICsgCiAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSBjb25mLmxvdywgeW1heCA9IGNvbmYuaGlnaCksIGFscGhhID0gLjI1KSArCiAgZ2VvbV9saW5lKGNvbG9yID0gJ2RvZGdlcmJsdWUnKSArIAogIGxhYnMoeCA9ICdJbmNvbWUnKQpnMiA9IAogIGdncHJlZGljdChtb2RfZ2FtMiwgdGVybXMgPSAiRWR1IikgJT4lIAogIGdncGxvdChhZXMoeCA9IHgsIHkgPSBwcmVkaWN0ZWQpKSArIAogIGdlb21fcmliYm9uKGFlcyh5bWluID0gY29uZi5sb3csIHltYXggPSBjb25mLmhpZ2gpLCBhbHBoYSA9IC4yNSkgKwogIGdlb21fbGluZShjb2xvciA9ICdkb2RnZXJibHVlJykgKyAKICBsYWJzKHggPSAnRWR1JykKZzMgPSAKICBnZ3ByZWRpY3QobW9kX2dhbTIsIHRlcm1zID0gIkhlYWx0aCIpICU+JSAKICBnZ3Bsb3QoYWVzKHggPSB4LCB5ID0gcHJlZGljdGVkKSkgKyAKICBnZW9tX3JpYmJvbihhZXMoeW1pbiA9IGNvbmYubG93LCB5bWF4ID0gY29uZi5oaWdoKSwgYWxwaGEgPSAuMjUpICsKICBnZW9tX2xpbmUoY29sb3IgPSAnZG9kZ2VyYmx1ZScpICsgCiAgbGFicyh4ID0gJ0hlYWx0aCcpCgooZzIgKyBnMyArIGcxICsgcGxvdF9sYXlvdXQobnJvdyA9IDIpKSAqIHRoZW1lX3RydWVNaW5pbWFsKCkKYGBgCgoKIyMjIDJEIFNtb290aHMKCmBgYHtyIG1vZF9nYW0zfQptb2RfZ2FtMyA8LSBnYW0oT3ZlcmFsbCB+IEhlYWx0aCArIHRlKEluY29tZSwgRWR1KSwgZGF0YT1waXNhKQpzdW1tYXJ5KG1vZF9nYW0zKQpgYGAKCgpgYGB7ciBtb2RfZ2FtM19wbG90fQojIHJlcXVpcmVzIHRoZSB2aXNpYmx5IHBhY2thZ2UsIG90aGVyd2lzZSBzZWUgdGhlIG1haW4gZG9jdW1lbnQ7IAojIGRldnRvb2xzOjppbnN0YWxsX2dpdGh1YigibS1jbGFyay92aXNpYmx5IikKdmlzaWJseTo6cGxvdF9nYW1fM2QobW9kZWwgPSBtb2RfZ2FtMywgbWFpbl92YXIgPSBJbmNvbWUsIHNlY29uZF92YXIgPSBFZHUsIHBhbGV0dGU9J2JpbGJhbycsIGRpcmVjdGlvbiA9IDEpCmBgYAoKIyMjIE1vZGVsIENvbXBhcmlzb24KCmBgYHtyIG1vZGVsX2NvbXBhcmlzb25fcmVkdXh9CkFJQyhtb2RfbG0yLCBtb2RfZ2FtMikKYGBgCgo=